knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
source(here('common_fxns.R'))
dir_str_mol <- '/home/shares/ohi/git-annex/impact_acceleration/stressors'
dir_str_10km <- file.path(dir_bd_anx, 'layers/stressors_10km')Collect stressor layers from CHI 2019 and reproject/aggregate them to the native resolution of the species range maps. This is a simple aggregation 11x as is done for the ocean, EEZ, MEOW, etc rasters.
Read in each raster, aggregate to a raster at the species projection and resolution. Because the projection is Mollweide for the source and target, and the base raster for species maps was created the same way, we can simply aggregate the original Mollweide CHI maps up by the same factor. For spp maps, the aggregation factor was 11\(\times\) to approximate 100 km2 cells: \(.934 \times 11 = 10.28 \text{ km}; \; (.934 \times 11)^2 = 105.6633 \text{ km}^2\).
Notes:
mazu:git-annex/Global/NCEAS-Pressures-Summaries_frazier2013/ice_mask_resampled (both .gri and .grd files needed)### List directories that end in "final"
str_dirs <- list.files(dir_str_mol, recursive = TRUE, full.names = TRUE) %>%
dirname() %>%
unique() %>%
.[str_detect(., 'final') & !str_detect(., 'archive')]
# [1] "archive/art_fish/final" "archive/comm_fish/final/dem_dest"
# [3] "archive/comm_fish/final/dem_nondest_hb" "archive/comm_fish/final/dem_nondest_lb"
# [5] "archive/comm_fish/final/pel_hb" "archive/comm_fish/final/pel_lb"
# [7] "archive/sst_old/final" "art_fish/final"
# [9] "benthic_structures/final" "comm_fish/final/dem_dest"
# [11] "comm_fish/final/dem_nondest_hb" "comm_fish/final/dem_nondest_lb"
# [13] "comm_fish/final/pel_hb" "comm_fish/final/pel_lb"
# [15] "direct_human/final" "land_based/final/nutrient"
# [17] "land_based/final/organic" "light_pollution/final"
# [19] "oa/final" "shipping/final"
# [21] "slr/final" "sst/final"
# [23] "uv/final" ### identify land-based stressors, and load coastal mask and sea-ice raster.
land_based_str <- get_str_cats() %>%
filter(category == 'land-based') %>%
.$stressor
coastal_str <- c('slr', land_based_str)
coastal_mask <- raster(file.path(dir_bd_anx, 'spatial/three_nm_offshore_mol.tif'))
seaice_mask <- raster(file.path(dir_bd_anx, 'spatial/sea_ice_mask_1km.tif'))
for(dir in str_dirs) {
### dir <- str_dirs[14] ### slr
### dir <- str_dirs[13] ### shipping
message('Processing stressors in ', dir)
files <- list.files(dir, pattern = '.tif$', full.names = TRUE)
# for(f in files) {
tmp <- parallel::mclapply(files, mc.cores = 12,
FUN = function(f) {
### f <- files[1]
str <- str_replace_all(basename(f), '_[0-9].+', '')
str_10km_file <- file.path(dir_str_10km,
str_replace(basename(f), 'mol', 'mol10km'))
if(!file.exists(str_10km_file)) {
message('Processing ', basename(f))
str_rast_orig <- raster(f)
### check for coastal or sst stressors
if(str %in% coastal_str) {
message('Masking ', basename(f), ' to coastal only')
str_rast_orig <- mask(str_rast_orig, coastal_mask)
}
if(str == 'sst') {
message('Masking ', basename(f), ' to remove sea ice')
str_rast_orig <- mask(str_rast_orig, seaice_mask)
}
str_rast_10km <- raster::aggregate(str_rast_orig,
fact = 11, fun = mean,
filename = str_10km_file,
overwrite = TRUE)
### check to make sure results match spp rasts
# ocean_area_rast <- raster(here('_spatial/ocean_area_mol.tif'))
# raster::compareRaster(str_rast_10km, ocean_area_rast,
# extent = TRUE, rowcol = TRUE, crs = TRUE,
# res = TRUE, orig = TRUE)
### TRUE!
} else {
# message('File exists: ', str_10km_file, '... skipping!')
}
})
}library(animation)
make_gifs <- function(map_stack, filename, layer_names = NULL) {
if(is.null(layer_names)) layer_names = names(map_stack)
x <- values(map_stack)[values(map_stack) > 0] %>%
min(na.rm = TRUE)
power_adj <- -1 * round(log10(mean(values(map_stack), na.rm = TRUE) * 2)) + 1
message('Adjusting by power of ', 1/power_adj)
map_stack_adj <- map_stack^(1 / power_adj)
ramp_breaks <- c(0, seq(x, 1, length.out = 100))
ramp_labels <- seq(0, 1, .2)
ramp_colors <- c('grey40', hcl.colors(n = 100, palette = 'viridis'))
capture.output({
saveGIF({
for(i in 1:nlayers(map_stack_adj)){
plot(map_stack_adj[[i]],
col = ramp_colors,
breaks = ramp_breaks,
axes = FALSE,
axis.args = list(at = ramp_labels^(1 / power_adj), labels = ramp_labels),
main = layer_names[i])
}},
interval = 0.5, movie.name = filename,
ani.width = 700, ani.height = 420)
})
return(invisible(NULL))
}str_10km_files <- list.files(file.path(dir_bd_anx, 'layers/stressors_10km'),
full.names = TRUE)
years <- 2003:2013
str_map_df <- data.frame(str_file = str_10km_files) %>%
mutate(str = str_replace(basename(str_file), '_[0-9].+', ''),
year = str_extract(basename(str_file), '[0-9]{4}') %>%
as.integer()) %>%
group_by(str) %>%
filter(min(year) <= 2003 & max(year) >= 2013) %>%
ungroup() %>%
filter(year %in% years)
stressors <- str_map_df$str %>% unique()
for(stressor in stressors) { ### stressor <- stressors[14]
### Animate the results
gif_file <- here('_setup', sprintf('figs/str_movie_%s.gif', stressor))
str_maps <- str_map_df %>%
filter(str == stressor) %>%
.$str_file
if(!file.exists(gif_file)) {
message('creating ', gif_file)
rast_files <- str_maps
map_stack <- stack(rast_files)
make_gifs(map_stack,
filename = gif_file,
layer_names = paste(stressor, years))
}
cat(sprintf('', gif_file))
}